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Abstract 

In this work we introduce the Object Kinetic Monte Garlo (OKMG) simulator MMonCa and simulate the defect evolution 
in three different materials. We start by explaining the theory of OKMG and showing some details of how such theory is 
implemented by creating generic structures and algorithms in the objects that we want to simulate. Then we successfully 
reproduce simulated results for defect evolution in iron, silicon and tungsten using our simulator and compare with 
available experimental data and similar simulations. The comparisons validate MMonCa showing that it is powerful and 
flexible enough to be customized and used to study the damage evolution of defects in a wide range of solid materials. 
Work submitted to Gomputer Physics Communications, 184 (12), 27032710 (2013). 
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1. Introduction 

The study of irradiation effects and defect diffusion in 
solid materials is a field of the maximum importance given 
its implication in technological solutions for the microelec¬ 
tronic companies and as structural materials for nuclear 
fusion and fission energy generation. Many materials have 
been studied under irradiation. For metals there are stud¬ 
ies on iron, tungsten [l| , cop^r-niobium , and others. 
For semiconductors silicon [^,0, silicon carbide 0 , ger¬ 
manium 0, gallium arsenide 0], and others. 

The physics involved within different crystalline solids 
when being irradiated is, to some degree, similar. Ini¬ 
tially, irradiation produces a population of simple point 
defects, typically correlated in interstitial vacancy pairs 
called Frenkel pairs. After some initial recombination of 
these pairs, one of the constituents of the pair might diffuse 
faster at certain temperature ranges. For instance, inter¬ 
stitials for iron around 130 K 0 or vacancies for silicon 
at room temperature llj. The moving particles agglom¬ 
erate around clusters that, at some point, might evolve 
into extended defects with different shapes and properties 
(l^ . In some cases, the extended defects are dislocation 
loops 01. Depending on the material, the extended de¬ 
fects might diffuse 141 or be immobile 0 ]. When diffus¬ 
ing, it is possible for them to react with other defects or to 
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reach the surface. When there are impurities present in the 
crystal, the impurities might diffuse in interstitial or sub¬ 
stitutional positions, and form clusters by agglomerating 
with other impurities of the same species, with intersti¬ 
tials and vacancies, or even with different impurities. In 
some cases, the role of such impurities is crucial to under¬ 
stand the behavior of the material under irradiation. For 
instance. He irradiation in metals, produces the formation 
of bubbles (He clusters) that are responsible for the change 


of the material mechanical properties in Fe [l6j, Gu [17 1 
W 01 and others. For semiconductor materials, cluster¬ 
ing of dopants is responsible for electrical de-activation of 
species like As 0 and B [^ . 

The study of such important phenomena through com¬ 
puter simulations has been a field of research for decades. 
First principle calculations are used to obtain the activa¬ 
tion energies and physical mechanisms of defect formation 
and diffusion [l9l [2l|. Lattice Kinetic Monte Carlo has 


been used to study macroscopic diffusivity |22|, cluster 
formation or recrystallization in heavily irradiated solids 


23|,[2^. Object Kinetic Monte Carlo (OKMG) is one of the 


preferred tools used to study defect evolution inside solids 
|25l427j| . And finally, when the internal micro-structure of 
defects is assumed not to be important and the concen¬ 
trations instead of the atomic positions offer enough infor¬ 
mation, finite element methods have been used. Many of 
these tools are well established and count with academic, 
open source, or commercial codes that are powerful and 
flexible enough to allow for fundamental research in all 
this wide range of materials. For instance, SIESTA ( 2 ^ . 
VASP 0 or Gaussian [0 are available for first prin- 
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ciple calculations. LAMMPS [3l|, GROMACS [S^ and 
many others are used to perform Molecular Dynamics, and 
there are many packages to run continuum (finite element 
method) simulations of which Abaqus and Ansys simula¬ 
tors, to name just a few, are well known and established. 

For OKMC, some existing codes are DADOS [1^, for 
diffusion of defects in silicon based materials, McDonalds 


33|, initially designed for silicon, Sentaurus Process KMC, 


a commercial software for Si based materials 3^, and 
LAKIMOCA [ 2 ^, a Lattice KMC used for simulation of ir¬ 
radiated metals. Nevertheless, there does not seem to be a 
clearly established, multi-material oriented, easy to access 
code for performing OKMC simulations. This lack could 
be negligible would it not have been for the extreme useful¬ 
ness played by OKMC simulations in the field of damage 
irradiation: being in the border between atomistic and 
continuum simulation, OKMC plays a very important role 
in using all the theoretical information on activation en¬ 
ergies obtained by the previously cited methods, and con¬ 


necting them to macroscopic experiments 2l|, [Sg . This is 


why in this work we want to introduce MMonCa, a recent 
OKMC simulator written in C-I--I- an integrated with the 
TCL script language, that wants to be multi-material, 
powerful, flexible and easy to use, filling the need for this 
type of codes that exist in the field of Monte Carlo simu¬ 
lation 371. 

This article is structured as follows: We will start re¬ 
viewing the KMC theory on Sec. [2] and the particular im¬ 
plementation of such theory on Sec. [3] Such implementa¬ 
tion will review the major modules of MMonCa: the time 
and space modules, (13.11 and 13.21 respectively 1 and the de¬ 
scription of all the implemented defect (object) types in 
Sec. 13.31 The results and validations are written in Sec. |4] 
starting with analytical calculations gTj and then iron 
(14.21) silicon (14.31) and tungsten (14.41) . Finally, we will sum¬ 
marize the work in Sec. [5l 


2. KMC theory 


The object Kinetic Monte Carlo algorithm goal is to 
follow the dynamic evolution of a system that might be 
out of equilibrium 38 - 4^ . It assumes that there are differ¬ 
ent states in the system, and that the transitions between 
these states are Markovian, that is, that the transition 
rates depend only on the initial i state and the final j 
state, and that such transitions are independent of time. 
These transitions are the input parameters of the al¬ 
gorithm. In our particular case, we model them assuming 
Harmonic Transition State Theory [4l| as Arrhenius laws 
with an activation barrier Eij (bigger than ksT for this 
approach to work) and a prefactor Pij: 


Tij — ^ij ^ ^^P( Eij /. 


( 1 ) 


The physical meaning of such barriers can be seen in 
Fig-in In such diagram Aij = Ej —E{ +Ak. The opposite 
Eji would be just the energy Eji = E^j. Assuming that the 



Figure 1: Energetic diagram for our KMC simulations, showing two 
states i and j and the formation and barrier energies related with 
them. 


concentration of particles in the i state is [i] and in the j 
state is [j], steady state will be reached when [i]rij = [j]rji. 
Using the notation stated in Fig.[T]and assuming Pij = Pji 
we have that reaching such state implies 




( 2 ) 


This relation does not include the barrier term Uk, 


but 

the difference in formation energies. Out of steady state 
OKMC can take care of the dynamic behavior of the sys¬ 
tem and provide a way to simulate time evolution. The 
inclusion of interacting particles is not always exact and 
implies some assumptions and limitations. Such assump¬ 
tions are, among others, a) a finite probability of trajec¬ 
tories to intersect without the reaction taking place for 
complex diffusion paths and/or complex object shapes, b) 
ternary reactions, c) collective movement and d) long-term 
reactions not happening or not being important, a) can 
be partly accounted for modifying the capture volumes, 
and b) and c) should not be a concern although there is 
a way to simulate collective movement in KMC ^ . Sim¬ 
ulating long-term reactions is possible through the inclu¬ 
sion of quasi-continuum fields, for instance Fermi-level for 
Coulombic effects in semiconductors l43ll a nd stress/strain 
computations for elastic interactions [44| . In this latter 
case, when the elastic interaction between particles is an 
important factor (for example in metals), it is usually in¬ 
cluded as a bias in the capture distance between different 
particles. 

Once all the transition rates for all the possible 
states in the system are known (that is, they are given as 
input parameters) the OKMC algorithm starts. For sim¬ 
plicity we will omit the initial state in the transition rates 
and write them as r^, being / the final state achievable 
from a particular initial state i. Using this notation, the 
KMC direct method [i^ is applied as follows: 
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1. Obtain the cumulative function 


i 


j=i 


(3) 


for z = 1, • • •, Being N the total number of tran¬ 
sitions in the given system. 

2. Compute two random numbers, r and s in the inter¬ 
val (0,1]. 

3. Find z, the event to perform, for which Ri-i < rR^ < 

R^■ 

4. Perform the event z: transform the particular chosen 
object from z to j. 

5. Increase the total simulated time by 


At = 


ln(l/s) 

Rn 


(4) 


6. Recalculate the affected rates. 

7. Return to step [T] until the requested physical time 
has been simulated. 


The above standard OKMC algorithm takes care of the 
time evolution only. Space dependence is intrinsic to the 
proper definition of each event. In our case the presence 
of physical defects that diffuse in space implies the need to 
include diffusion as a transition rate, and to define algo¬ 
rithms for space migration and particle interaction. Con¬ 
sequently, our OKMC simulator for damage evolution in 
solids contains the following modules: 

• Objects (defects) and the list of their associated tran¬ 
sition rates and actions. 

• A rate manager to compute time evolution and to 
pick up the event to perform. 

• A space manager to manipulate space translations, 
neighbor search and defect interactions. 


3. Implementation 

Fig. [5] shows the overall structure of our simulator. 
MMonCa has been implemented as C-I--I- extensions of the 
TCL [ 3 ^ language. This allows us to use an already ex¬ 
isting and well known language for the input script and 
to implement a user interface. The commands that have 
been extended allow the user to define a simulation cell, 
the 3D definition of the material structure of the simula¬ 
tion, reading of damage from an external file, annealing 
the damage, reading and writing the parameters needed 
for the simulation from the input file, and output of differ¬ 
ent quantities generated during the simulation, being the 
most important the concentrations and defect position and 
types. The rest of modules are described below. 



Figure 2: Overall structure of the MMonCa simulator. The user in¬ 
terface relies on a layer of the TCL interpreter, extended to support 
OKMC. The extension relies on specialized modules to control the 
space, time and defects. Several defects are supported as the ob¬ 
jects to be simulated: extended defects (ED), mobile particles (MP), 
damage clusters (DC), multi clusters (MC) and interfaces (Int). 


3.1. Computation of time: rate manager 

Fig. |3] shows graphically the idea behind the event se¬ 
lection involved in step [3] of the OKMC algorithm previ¬ 
ously explained. Once an updated list of all the transitions 
associated with the objects being simulated is generated, 
one of them is chosen proportionally to such rates. The 
At = — ln(s)/i?Ar associated with the simulation of such 
event is independent on the event chosen, depending only 
on the whole system. In practice, iterating through all the 
cumulative rates to find the one to be performed is not effi¬ 
cient when there is a large number of rates. For this reason, 
our simulator does not contain a transition bar with all the 
rates, but rather a binary tree, where the access time to 
each rate is not proportional to the number of them N 
but to log 2 (V). In the one hand, this improves the access 
time to the chosen event, on the other hand, the binary 
tree degrades the insertion, deletion and modification time 
for rate insertion from a constant time to also oc log 2 (V) 
time. Overall, the balance is positive when there are a 
large number of rates in the system and more selection of 
rates than insertions, modifications or deletions. 

3.2. Space organization and neighbor location 

The space is divided in small prismatic elements us¬ 
ing a tensor mesh. Space is assumed to be homogeneous 
(material, temperature and other fields) inside each small 
element. Each element obtains its material definition by 
calling a user-defined procedure that allows the specifica¬ 
tion of the material structure in the simulation. This way, 
very complex shapes containing different materials can be 
simulated. When two consecutive mesh elements have dif¬ 
ferent materials an interface object, as shown in Fig. 0^), 
is built between them. 
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Figure 3: The OKMC algorithm contains a list of all the transitions 
associated with the objects being simulated, and picks the next one 
proportionally to such rates. That can be seen graphically as getting 
a random number uniformly distributed in [0,Rn) s,nd picking up 
the event “aligned” with such number. 



Element 


Face 


Figure 4: Space is divided into small prismatic elements (rectangles 
in this 2D representation) using a tensor mesh, a) An interface is the 
union of all element faces between adjacent different materials, b) 
The capture distance Tc of every particle is defined independently, c) 
The capture distance of clusters is built as the union of the capture 
distances of their constituent particles. 


Efficient neighbor search is implemented by having a 
standard link cell method [i^ . Once the list of neighbors is 
obtained a look-up table is used to implement user-defined 
allowed interactions. 

The capture distance he, shown in Fig. 0 Jd), must be 
provided for every single particle. It is typically of the 
same order as A, the microscopic migration distance. In 
our simulator, any non point defect (except interfaces) 
is created by the agglomeration and tracking of its con¬ 
stituent particles. This implies naturally that the cap¬ 
ture distance of any defect is the overlay of all the capture 
distances of all the defect constituent particles as seen in 
Fig. lib). It also means that extended defects can have any 
shape and a capture volume that will adapt to it as long 
as the particles are configured to form such shape. 


can be introduced on purpose to improve the material fea¬ 
tures (for instance, doping of semiconductors to produce 
devices). In any case, we define the objects of our simula¬ 
tor as the defects introduced in the material. In particular, 
we classify such defects as interfaces (Int), mobile particles 
(MP), damage clusters (DC), extended defects (ED) and 
multi-clusters (MC). The properties we simulate for each 
of them are described next. All these are considered ob¬ 
jects of the OKMC simulator and are treated in a similar 
way. Each object, to be included in the simulator, needs 
to have the following data and functions defined: 

• Number of events associated with the object. For 

instance, three for MPs: diffusion, breaking-up and 
creation of Frenkel pairs {A ^ A +1+ V —>■ -|- 

V{I)) to react with impurity atoms. 

• Rate associated to each event. In our example for 
MPs, computing the diffusion, breaking-up and in¬ 
jection rates, or returning zero if they do not apply. 

• Functions to perform each event when it is chosen by 
the OKMC algorithm. In the example, an MP needs 
the implementation to move the particle, break it up 
or create and/or react with Frenkel pairs. 

Some of the explained events (break up and creation/reaction 
with Frenkel pairs) implement reactions similar to AB —>■ 

A + B (for instance, Ci ^ C +1 oi Hcg —>• Hci + V). The 
forward reaction A + B AB is implemented through 
diffusion. For this forward reaction to happen two things 
are needed: a) A moving to the neighborhood of B, or 
B into A, and b) the reaction being allowed. Diffusion is 
implemented as an event for all defects but interfaces. At 
the end of such event, a look for neighbors is performed to 
detect potential reacting species as explained in Sec. 13.21 
To properly react with such species, two more algorithms 
are needed in each KMC object 

• A look-up table that establishes whether the reac¬ 
tion is possible or not (taking into account possible 
reaction barriers) 

• A function that implements the interaction itself, 
taking the reactants and transforming them in the 
result. 

Finally, since during reactions the reactants are de¬ 
stroyed and the result is created, each object requires a 
constructor and a destructor that is able to properly build 
and erase respectively the objects from the KMC simula¬ 
tor. 


3.3. Defect structure 

We want to apply the theory of OKMC to the partic¬ 
ular problem of simulating the evolution of damage intro¬ 
duced into a solid. Such damage can be introduced as an 
undesired side effect of the material application (for in¬ 
stance, when using the material in a fusion reactor) or it 


3.3.1. Int: Interface 

Fig.|T]shows how the plane between two different mate¬ 
rials or a material and the outside world is defined using an 
interfaces object. Interfaces can create and inject MPs (Is 
and Vs, and also impurities that were previously trapped). 
These emissions can be done to either side, assuming the 
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Figure 5: Three phase segregation model. Particles can be captured 
and emitted at either side, but the binding energies, migration ener¬ 
gies and capture barriers might be different at each side. 


z/Q where the activation energy for 

such example would be set as Ef{I)+Ef{V)—Ei,{Hei)+ 
Em{I), being Ef the formation energy. 

MPs can interact with each other to form more complex 
defect objects: for instance / + / producing DCs or EDs, 
or HeV + He producing MCs. 


3.3.3. DC: Damage cluster 

DCs are irregular agglomerations of I and V with a 
non-instantaneous recombination rate. This mechanism 
simulates the recombination time needed by IV pairs in 
some systems, that although small, is not null, to anni¬ 
hilate both defects T7|. The rates associated with DCs 
are: 


MP may exist there. In the particular case of impurity 
emission, the model implemented corresponds to a three 
phase segregation model. Such model is shown in Fig. [S] 
MP impurities can be at the interface by overcoming the 
barrier to reach the interface (Ebarrier + Em). Then, they 
have a rate v = r'g exp(—Eemit(side)// cbT) to be emitted 
to either side. Eemit is set as Ef,-f Ebarrier + ■ Similarly 

to equation [21 it is easy to see that the segregation coef¬ 
ficient, defined as the ratio between the concentration of 
particles at both sides at equilibrium, is 

S = eMiEb^) - Eb{2))/kBT). 

When any diffusing defect arrives at the interface it 
can be annihilated according to certain probability set by 
the user. This applies to MPs, EDs, DCs and MCs. 


1. Recombination of IV pair with 

ly = Vo ex.p{—Eiv(size)/ksT). 

2. Emission of MPs. The constituent particles can be 
emitted with a rate 

V = i/oexp(-£'emit(size)/fcB7’) 

until the cluster dissolves. The activation energy for 
emission Eemit (size) is computed as the binding en¬ 
ergy for each size plus the migration of the emitted 
particle. 

3. Transformation into an ED. The transformation rate 
is computed as 

Vq exp( Etransforin(size)//c^T). 

4. Diffusion by random walk with rate 

V = z/oexp(-E,.„(size)/fcBT). 


3.3.2. MP: Mobile particle 

Single [He, /, V) or paired defects (CV, Cl) are de¬ 
fined as MP in our simulator, where paired interstitial de¬ 
fects are assumed to be the same as impurities in the in¬ 
terstitial position [Ci = Cl). The transitions associated 
with these MP objects are: 

1. Migration, by simulating the random walk of small 
diffusion events with fixed migration distance A in 
one of the three perpendicular axes of the system, 
randomly chosen for each jump. The migration rate 
for mobile particles is computed asv = vo exp){—EmlkB 
where vq and Em are the input parameters for mi¬ 
croscopic diffusivity. 

2. Break up of a pair (or kick off mechanism) of I or 
V impurities. For instance, CV C + V. The 
break-up frequency equals v = z/q exp(—E^fc/fesT) 
with Ebk being the activation energy for break-up. 
Such activation is computed as binding energy plus 
migration energy of the emitted particle. 

3. Injection of extra Is or Vs by creating an IV pair, 
capturing the I or the V and emitting the other 
(also called Franck-Turnbull mechanism). This re¬ 
action applies for instance to He substitutional in 
W: Hcs ^ Hei + V. Its rate is modeled as v = 


Rates 12112] and m are non-null when the damage cluster 
contains only Is or Vs, but not both. DCs can also interact 
with MPs and EDs. 

3 . 3 . 4 . ED: Extended defect 

EDs are agglomeration of interstitials (/„) or vacan¬ 
cies (Vn) with particular shapes that can emit their con¬ 
stituent particles, transform into other EDs, migrate and 
trap/detrap impurities that might stop their diffusion. In 
contrast with DCs, EDs contain only Is or Vs but never 
j^both. EDs can adopt different shapes to adapt to the real¬ 
istic morphology of extended defects in different materials. 
In particular, they can be defined as a) planes (similar to 
{311} defects in Si h ) disks (similar to dislocation 

loops in Fe 21| or Si [^), c) spheres (voids in Si and 


other materials 0 ) and d) irregular clusters (no special 
shape). 

The transition rates defined for the different events are: 

1. V = i/Q exp(—Eemit(size)/fcBr) for emission of MPs, 
being Eemit (size) the addition of binding energy plus 
migration energy of the emitted particle. 

2. V = z/Q exp(— i?transform(size)/fcBT) for transforma¬ 
tion into other EDs, being each transformation acti¬ 
vation energy and prefactor defined by the user. 
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i. V = vq /ksT) for migration. 

A. V = i/Q exp(—i?detrap(particle)/fcsT), one value for 
all sizes, to detrap previously captured particles. In¬ 
dependently the user can specify whether the trapped 
particles are only decorating the extended defects or 
also stops its diffusion. In this latter case, the detrap¬ 
ping rate plays a very important role for the overall 
diffusion rate of extended defects in the presence of 
traps. 

EDs can react with new incoming MPs. Depending on 
the nature of the incoming particle, the defects will grow 
(/„ -!-/—>■ /„+i), annihilate it instantaneously (/„ -I- E — > 
or trap it (/„ -|- C —>■ CIn)- Reactions with other 
extended defects are also permitted {In + dm —>■ In+m, 
In + Vm, Vm-n assuming m > n). They can also react 
with MCs (assuming that the final product is defined) and 
DCs and transform into a different ED with the same size 
(for instance, /„ < 111 >—>■ < 100 >). 


3.3.5. MC: Multi-cluster 

MCs are the agglomeration of several impurities with 
either Is or Vs. They can play different roles in the physical 
systems under consideration. They allow the simulation of 
helium cluster formation in metals like Fe [H^ and W [5l| . 
In semiconductors, clusters of dopants with interstitials 
and vacancies deactivate partially the implanted dopants 


by forming agglomerations like As 4 V 
interstitial clusters 


52, 


or boron 


The different events that MCs can perform are: 
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Figure 6: (Color online). Comparison between KMC (symbols) and 
the theoretical solution (lines) for the time evolution of a spatial 
concentration of non-interacting particles with periodic boundary 
conditions. 


instance, two rates are needed for He 2 V —>■ He HeV, 
once considering the barrier of Em{He) and the other 
En^iHeV). 

MCs can react with MPs, EDs and other MCs as long 
as the formation energy of the final result is included as a 
parameter. Even in those cases, a probability to reject the 
reaction 


1. Emission of their constituent particles as MPs. The 
activation energy for emission is computed as the 
formation energy difference between the final and the 
initial state when positive, plus the migration energy 
of the emitted particle. The potential energies of all 
the clusters are required input parameters for the 
simulation. 

2. Emission of constituent particles in pairs. For in¬ 
stance, AnVm An-iVm-i + AV. The activation 
energy being Ef{An-iVm-i) + Ef{AV) -Ef{AnVm) 
when positive, plus Em{AV). 

3. Injection of non-existing Is or Vs by Frenkel pair cre¬ 
ation {Anljn ^ AnIm +1 -f V Or AnMrn ^ AnVni+l^- 
The activation energy equals to Ef{AnIm +i)+Ef{I) + 
Ef{V) — Ef{AnIm) when such value is positive (zero 
otherwise), plus Em{I) for injection of Is. 

4. Migration. The migration rates are defined in a sim¬ 
ilar way to all the other migrations as 

vq exp(—i?m(cluster)/fcBT). 

The number of different MP emission mechanisms for 
a simple MC can be high. For instance, He 4 V clusters can 
make transitions to He 4 -I- V, He^ -I- HeV, He^V He, 
He 4 V 2 1 and + -ffel. For clusters breaking into 

elemental particles (MPs) the simulator also has to con¬ 
sidered the migration energies of both constituents. For 


P = exp[(E}-E/)/fcBT] 

is defined to account for the barriers involved in the for¬ 
mation of the new cluster. If Ej < Ej the reaction always 
happens. 

4. Results 

This section describes the validation of the code by 
comparing with theoretical values and experimental re¬ 
sults or other simulations in three different materials: iron, 
silicon and tungsten. 

4-.1. Theoretical results 

Fig.lHshows the comparison between KMC simulations 
(symbols) of the temporal evolution of an initial distribu¬ 
tion of particles and the exact, theoretical results (lines). 
The calculation of the theoretical results has been done 
similarly to Ref. [H^. The KMC simulations have been 
run for 10 seconds in a 100 x 300 x 300 nm^ simulation cell 
with a total number of 448820 particles. The diffusivity of 
each particle was set to lOOnm^s”^. Further comparisons 
with theoretical results, not shown here, have been done 
for interacting particles (for instance, diffusion of impuri¬ 
ties through intermediate species A 1 Ai and break 
up), reaction with interfaces or sinks, correct establishing 
of equilibrium concentrations, etc. 
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Table 1: OKMC Iron model 
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Figure 7: (Color online). MMonCa simulation of the evolution of de¬ 
fects and resistivity recovery during isochronal annealing. Top figure 
shows the derivative of the total defect concentration (red curve) 
being compared with experimental results M for recovery stages 
(black arrows). Bottom figure shows the total simulated concen¬ 
tration of defects and the different defect contributions (lines) during 
the isochronal annealing of the sample after electron irradiation com¬ 
pared with previous theoretical work (symbols) |21| . 


4-2. Iron 

The study of defect kinetics in irradiated iron is a prob¬ 
lem of primary importance for the aging of materials in 
the nuclear industry. Experimental work has been done 
by resistivity recovery experiments in high-purity electron- 
irradiated iron by Ref. [l0|, with irradiation doses in the 
range ft; 2 x 10“® to ~ 200 x 10“® displacements per atom 
(dpa). In these experiments, the resistivity of the metal is 
recorded during an isochronal annealing. The derivative 
of the resistivity versus the temperature shows clear peaks 
that are called recovery stages. These stages are related 
to different physical mechanisms involving the recombi¬ 
nation, migration, growth and dissociation of the defects 
formed during irradiation and subsequent annealing. In 
particular, five important stages have been detected for 
iron. 

• Stage Id 2 , observed at 107.5K related to the recom¬ 


Object 

Migration 

Species 

Parameters 

MB 

Yes 

/ and V 

Ref. [21, Ml 

ED 

Yes 

In small clusters 

Ref. [21, Ml 

ED 

Yes 

<111> In clusters 

Ref. 

ED 

size < 5 

Vn clusters 

Ref. [^,M] 


bination of Frenkel pairs. 

• Stage Ie around 123 to 144 K, as the result of the re¬ 
combination of I and V belonging to different Frenkel 
pairs through the migration of interstitials. 

• Stage II is suggested to happen when the I 2 starts 
to diffuse, around 164 to 185 K. 

• Stage III attributed to migration of Vs, around 220 
to 278 K. 


• Finally stage IV, around 520 to 550 K produced by 
the dissociation of defect clusters formed during the 
previous stage III. 


Fig. [3 shows the simulated isochronal annealing of 2 x 
10 “^ dpa irradiated iron, together with the experimental 
stages (black arrows). It can be seen that the agreement 
with experiments [l0| and with previous simulations done 
by other groups is good dllEi, especially taken into ac¬ 
count that the compared results are produced by two dif¬ 
ferent KMC methods (Event versus Object). A brief sum¬ 
mary of the models and parameters used for such simula¬ 
tion is shown in Table. [T] 


4-3. Silicon 

The evolution of defects in silicon has been a subject 
of intense research for the past decades. Its interest relies 
on the need of semiconductor manufacturers to understand 
the Si system to produce more powerful electronic devices. 
One particular subject of study has been the characteri¬ 
zation of damage by Si implantation. The evolution of 
such system contains many phases that are nowadays well 
known . The initial implantation produces a high pop¬ 
ulation of Is and Vs, where the V diffuses, even at room 
temperature implantations. During this initial stage Is and 
Vs do not recombine instantaneously, and tend to form 
DCs of various sizes. Depending on the particular implan¬ 
tation conditions, the amorphous pocket population might 
in some cases grow big enough to partially amorphize the 
sample. In other cases, dynamic annealing of the gener¬ 
ated damage, that is, the annihilation of IV pairs during 
a cascade and the next one, might be enough to avoid 
amorphization. 

Once the implantation has finished, the sample is pro¬ 
cessed to anneal out the defects. This typically eliminates 
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Figure 8: Interstitial supersaturation as a function of time after a 
40keV, 2 X 10^® cm“^ Si into Si irradiation at different tempera¬ 
tures 600, 700 and 800 °C. Symbols: Experimental data taken from 
Ref. [T^ . lines: simulation results using the OKMC code MMonCa 
presented in this work. 


Table 2: OKMC Silicon model 


Object 

Migration 

Species 

Parameters 

Int 


I, V creation 

Ref. [M] 

DC 

No 

InVn 

Ref. [34, 47| 

MP 

Yes 

I and V 

Ref. [57, Ml 

ED 

No 

{311} In clusters 

Ref. 

ED 

No 

Vn voids 

Ref. [^ 


all the DCs, leaving only small extended defects in the 
beginning. Such extended defects are composed of the 
extra interstitials introduced by the implantation. During 
the annealing, the small, irregular interstitial clusters emit 
their constituent particles. This produces an almost con¬ 
servative Ostwald ripening where big defects grow at the 
expense of small ones. At some point, the defects are big 
enough to be seen through the microscope, getting a char¬ 
acteristic {311} shape. Further annealing of these defects 
produces its dissolution or the formation of the very stable 
dislocation loops. 

Figure [5] represents the comparison of experimental 
supersaturation (concentration of Is in equilibrium ver¬ 
sus measured concentration) with the simulated results of 
MMonCa. The experimental results are taken from Ref. 121. 
In the experiment an implantation of 40 keV, 2 x 10^^ cm“^ 
Si+ into Si was followed by annealing at 600, 700 and 
800 °C. Table [5] shows the objects that we have defined 
and the references we use for the correct parametrization 
of such objects. Excellent agreement with both experi¬ 
mental data 121 and simulations is achieved. 


Table 3: OKMC Tungsten model 


Object 

Migration 

Species 

Parameters 

MB 

Yes 

I and V 

Ref. 

MB 

Yes 

He 

Ref. 

ED 

Yes 

In 

Ref. [^ 

ED 

Yes 

Vn 

Ref. [^ 

MC 

Yes 

HCn 

Ref. [^ 

MC 

No 

HCnV^rn 

Ref. 

MC 

No 

HCnIn 

Ref. 

MC 

No 

CVn (traps) 

Ref. 

MC 

No 

CIn (traps) 

Ref. [^ 


4-4- Tungsten 

Tungsten is usually proposed as an appropriate mate¬ 
rial for nuclear fusion reactors due to a number of features: 
low-activation, high melting point, low sputtering yield, 
high thermal conductivity and low thermal expansion. W 
is proposed as armor material for inertial confinement fu¬ 
sion by laser with direct drive targets [H^ . For future mag¬ 
netic fusion power plants W is considered the material of 
choice for the first wall and divertor [^. Consequently, 
simulation of irradiation-induced damage in W by OKMC 
can help in the understanding of such a material 611. 

The parametrization used to model W has been taken 
from Ref. [5l| and is summarized on Table [S] It consti¬ 
tutes a complex model that lets all defects interact with 
each other and with traps and allows for cluster formation. 
All pure clusters may migrate. In the particular case of 
interstitial clusters the migration is ID along <111> di¬ 
rections. Simulation boxes of dimensions 399 x 400 x 1001 
in lattice units, with lattice parameter A = 0.317 nm were 
used. The boundary conditions were periodic for y and z. 
The X surfaces (both) were assumed to allow the desorp¬ 
tion of incoming defects with a probability of 100%: all 
approaching defects are annihilated. 

We compare our results with those of Becquart and 
co-workers for the amorphous case 6^, i- e., we ignored 
the crystal structure of W when calculating the Frenkel 
pairs created by every incoming ion. 100 appm of C were 
introduced as static traps acting on interstitials and va¬ 
cancies, as well as on their clusters. We used the same 
irradiation conditions as Becquart (3keV He irradiation 
at 5 K and 16 He per second up to a dose of 12 ppm). 
We realized during the validation of MMonCa that the re¬ 
sults strongly depend on the initial conditions {He, V and 
/ distributions). Therefore, we used the same initial de¬ 
fect distributions as Becquart and co-workers obtained for 
amorphous W (^ . After the implantation the tempera¬ 
ture was decreased to 1 K and isochronal annealing steps 
of 2 K for 60 s were simulated. 

Fig- m compares the results of MMonCa (lines) with those 
presented in Ref. (symbols) concerning the evolution 
of interstitials, vacancies and helium remaining in the sim¬ 
ulation box. Fig. [TUI on the other hand, displays the num- 














ber of trapped and free helium atoms remaining after every 
annealing step. We can observe that the agreement is fair 
over the whole simulation for the different types of defects. 
Also defect clustering (not shown) is fairly reproduced. 
However, some discrepancies appear. We mainly attribute 
them to the different procedures employed by the codes 
to account for defect trapping. The code used by Bec- 
quart and co-workers considers that: (i) every defect has 
an associated capture distance; (ii) clusters are spherical 
objects with an associated capture distance that in general 
grows with the number of constituents; and (iii) whenever 
the capture volumes (defined by the capture distance) of 
two defects overlap, trapping occurs. On the other hand, 
MMonCa associates a capture distance to every single defect, 
whereas the clusters are formed by the agglomeration of 
single defects in different configurations (see Section [33) • 
In any case, the identity of the single defects is not lost 
and trapping occurs when an object falls within a distance 
smaller than the capture distance of any single defect. 
Therefore, the trapping procedures are different and this 
turns out to be the source of the small discrepancies found 
when comparing the results. Note that the capture dis¬ 
tances used in Table 5 of Ref. [sij can not be directly used 
in MMonCa because trapping is defined in different ways. 
In principle we must use capture distances approximately 
twice the size than those previously reported by Becquart 
and co-workers to account for their trapping criterion. We 
have found that the best results are obtained when we mul¬ 
tiply the capture distances given in Table 5 of Ref. ^ by 
2.3 for the mobile particles /, V, Hel and HeV and by 1.5 
for He, C, Cl and CV. With these values, MMonCa slightly 
overestimates the interstitial loss at 7 K and the helium 
release at around 300 K Fig. (jO]). In addition, the helium 
trapped fraction at low temperatures (Fig. (TO]) turns out 
overestimated (the helium free fraction is underestimated). 
The different trapping procedures used in both codes are 
responsible for slightly different cluster formation during 
implantation. This, in turn, has consequences for the final 
evolution of the defects during the isochronal annealing. 
The largest differences between the codes are related to 
the evolution of big clusters, because the optimization of 
the capture distances can not account for the values as¬ 
signed to every cluster by Becquart. However, despite the 
small discrepancies observed, we conclude that MMonCa is 
able to reproduce complex results according to the expec¬ 
tations. 


5. Conclusions 

In this work we have reviewed the simulation tech¬ 
niques of the evolution of damage in irradiated solids and 
we have introduced the OKMC simulator MMonCa and ap¬ 
plied it to show the defect evolution in three different ma¬ 
terials. We have started by explaining the theory of KMC 
and showing some details of how such theory has been im¬ 
plemented by creating generic structures and algorithms in 



Figure 9: (Color online). Comparison of the number of interstitial, 
vacancy and helium as simulated in this work (lines) and in Ref. |62|| 
(symbols). 



Figure 10: (Color online). Comparison of the number of total, 
free, and trapped He atoms as simulated in this work (lines) and 
in Ref. (symbols). 
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the objects that we want to simulate. We have then repro¬ 
duced experimental and simulated results in iron, silicon 
and tungsten using our simulator. The different compar¬ 
isons show that MMonCa can be successfully used to study 
the damage evolution of defects in solid materials validat¬ 
ing the OKMC approach and the particular implementa¬ 
tion into the MMonCa simulator, that we hope will be of 
help for the materials research scientific community. 

A copy of the simulator described in this work can be 
obtained at the following web page: 

http;//www.materials.imdea.org/MMonCa/ 
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